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We have set up and tested a pipeline for processing the data from a spherical 



O 

■ gravitational wave detector with six transducers. The algorithm exploits the multi- 



channel capability of the system and provides a list of candidate events with their 
arrival direction. The analysis starts with the conversion of the six detector outputs 



(N 

^NJ ' into the scalar and the five quadrupolar modes of the sphere, which are proportional 



to the corresponding gravitational wave spherical components. Event triggers are 
then generated by an adaptation of the WaveBurst algorithm. Event validation and 
direction reconstruction are made by cross-checking two methods of different inspira- 



CN ' tion: geometrical (lowest eigenvalue) and probabilistic (maximum likelihood). The 

, combination of the two methods is able to keep substantially unaltered the efficiency 

I and can reduce drastically the detections of fake events (to less than ten per cent). 

lO ' We show a quantitative test of these ideas by simulating the operation of the reso- 

o ; 

^0 . nant spherical detector miniGRAIL, whose planned sensitivity in its frequency band 



(few hundred Hertz's around 3kHz) is comparable with the present LIGO one. 



^ ! PACS numbers: 04.80.Nn,95.55.Ym 



I. INTRODUCTION 



What makes a resonant spherical detector 



2j really different from bars [3|] and inter 



ferometers Q] is the fact of being a multichannel detector. As result, a sphere has almost 
isotropic sensitivity and enables to reconstruct the gravitational wave (GW) direction. After 
the seminal work of Wagoner and Paik JsJ, where the basic features of the detector have been 
pointed out, more and more detailed sphere models and configurations have been studied by 
several authors 0, Q, [s, 9, 10, 11, 1^, Is], M\- In particular, various solutions have been pro- 
posed to the problems of the ideal transducers configuration [12| , parameters reconstruction 



0,1111], and multidimensional data analysis [9|. 

The analysis method that we have set up and studied in the present work is the result of 
a synthesis and a refinement of some of these contributions, and is intended to represent the 
core of the pipeline that we are building for the miniGRAIL detector. MiniGRAIL is one of 
the two small prototypes of resonant spherical detectors which will be soon operating in the 
3kHz reg ion, with the possibility of detecting quasi-normal modes of strange core neutron 



stars 



15l | or burst from sub-solar mass black hole coalescences 16|. 



A greater variety of sources could be studied by building a larger sphere operating in 
the sub-kHz region and our method, being completely general, could be applied to such 
detectors as well. 

Sec. HI] of the paper is devoted to building an accurate detector model and to the gener- 
ation of a simulated set of data. We concentrate on the case of six transducers placed in a 
specific configuration (the so-called TIGA [10]), which best preserves the spherical symme- 
try of the detector, and has been adopted by the miniGRAIL team. 

For the same reason as above, the parameters of the numerical sphere model (mass, radius, 
Poisson's ratio, modes and transducers frequencies, temperature, quality factors, readout 
electronic components) have been chosen in order to mimic the features of miniGRAIL, 
including asymmetries and possible "imperfections" . 

We work out the spectrum of each transducer's output and we generate simulated data 
starting from the elementary noises which are involved at different stages in the detector. 
Then, we dicuss the problem of extracting the GW parameters out of the data. 

In sec, mil we describe the method to generate trigger of events. We used WaveBurst, the 
event trigger generator used by the Burst working group of the LIGO- VIRGO joint collabo- 
ration, suitably adapted for our multi-mode analysis. Once an event trigger is obtained, the 
five quadrupolar modes provide a redundant description of a GW signal, which depends on 
four parameters (the amplitudes of the two polarizations and the two angles identifying the 
source direction). We exploit this redundancy to distinguish true GW signals from other 
non-Gaussian noise (like glitches, for instance). The multimode analysis is then able to 
compute all the relevant quantities of a GW and to reduce the false alarm rate. 

Sec. IIVI contains some results of the analysis, showing the level of signal-to-noise ratio 
required to obtain the desired accuracy in the direction reconstruction and efficency. 
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Figure 1: The readout scheme of one of the tranducers. The surface of the sphere is modeled 
as an oscillating mass M, whose displacement is amplified by the six trasducers of mass m^, 
k = 1..6, placed at different positions on the sphere. For each transducer an electric circuit 
(primary) converts the mechanical displacement into a current Ik which is then converted into a 
larger current Ik in the secondary circuit. A SQUID finally amplifies the secondary current. 

II. MODELING THE DETECTOR 



A. Equations for modes, transducers and readout currents 

The sphere is well modeled [lo|, Q as a set of coupled mechanical and electric oscillators, 
describing the dynamics of the relevant sphere vibrational modes, of the transducers, and 
of the electrical circuit that are at the core of the readout devices, see fig. [TJ 
As the relevant equations have been already discussed in detail by several authors (see for 
instance the references above), we jump to the mathematical core of the problem, skipping 
introductory material and definitions that can be found in the literature. 

In Fourier space, the detector is described by a system of linear and algebraic equations 



{^N, Qk, Jk, h} — Z /tv, fl, fk, fk, fk 



(1) 



where the are the amplitudes of the various radial modes of the sphere, the Ik, Jk are the 
currents flowing in the k—th transducer electric circuit of the readout, the qk the positions 
of the relative mechanical oscillators, and the /'s the stochastic forces related to the various 
dissipative components of the detector. As to the operator Z, its exact form can be deduced 
from the equations contained in 10|, llJ] : for our purposes, what matters is that this operator 
can be computed numerically as a function of the detector parameters. 
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Figure 2: The simulated output of transducer ^^0. 



Since the noise spectral densities of the stochastic forces /'s are known (again, see 



Id, 



14|). a possible detector configuration can be produced thanks to a random number 



generator; then the corresponding output currents spectral densities can be derived through 
eq. (PP) , as shown in figure [2J 

B. Prom currents to the GW wave modes 



When a GW impinges the detector, a deterministic force is added to the quadrupolar 
components of fjy: 



/jv ^ /jv + ^RxNU}^hN for 



7V=0,4 



0, 



A. 



(2) 



(3) 



where R is the sphere radius and = ymN,ij^^^'' ^ being ra* the versor of the arrival 
direction of the wave and rriN = 0, Ic, Is, 2c, 2s, according to the conventions used in j^. 

When six transducers are present, as in the TIG A configuration, the system is in principle 
overconstrained as the 5 quadrupolar Hm's are to be determined by the 6 outputs. 

In the idealized case, the solution is provided by the mode channels , i.e. 5 linear combi- 
nations of the currents which are in one to one correspondence with the sphere quadrupolar 
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Figure 3: The mode /iq as it has been reconstructed starting from an actual noise realization, 
superimposed with its expectation value. 



modes lOj, while the sixth, remaining combination is insensitive to any such modes. 

In the presence of noise and of asymmetries however the exact definition of mode channels 
is not trivial and require the addition of a sixth vibrational mode of the sphere in the problem, 
in order to make the system invertible 



We chose then to include the scalar mode in the system, as this is the mode that in the 
idealized case is probed by the sixth linear combination of current outputs. Moreover, the 
determination of the scalar mode allows, in principle, to test alternative theories of gravity 
(where scalar GW's can exist), or to build a veto, as events exciting such mode cannot be 
reduced to Einstein general relativity GW's. 

Given these premises, the problem can now be numerically solved (for every value of u) 
to give 



(4) 



where T^^i^) defines the transfer function (which in this case is actually a transfer matrix) 
of the system. 

We can now use eq. (jlj) to express any realization of the detector output in terms of GW 
quadrupolar modes, as shown in figure [3l Moreover eq.(jl]) can be combined with the results 
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of the previous paragraph in order to give the spectral density matrix: 

< 4(cj)/*,(cj') >^ Tjvfc ^< hNiuj)h%,{iu') >= S'^N,{iu)5{iu - lu') , (5) 
which can be used to estimate the detector sensitivity. 

III. THE ANALYSIS PIPELINE 

We have just hinted how to generate six time series corresponding to the quadrupolar plus 
scalar mode. As at this stage we are not interested in the scalar GW component, then we 
will restrict the analysis to the five quadrupolar modes, which will be enough and necessary 
to reconstruct the most general symmetric, traceless 3x3 matrix. 

Still this is a redundant description of a GW, and we will exploit this redundancy to 
discriminate between real GW signals and excitations of the modes due to disturbances 
other than gravitational. 



A. The scalar trigger 

If one knew exactly the form of the expected signal, the optimal strategy would be to 
perform a multidimensional matched filtering, but we would like to find triggers out of the 
stretch of data, without any need of a detailed knowledge of the signal. This can be done 
by building the following quantity: 

H{u;) = {h-h)s, (6) 

which is a generalization of the one proposed in [L^, and takes account of the differences 
and of the correlations among the various quadrupolar modes, as the r.h.s. is defined by 

{X-Y)s^ E xUu^){SU,)-\u)Yr,iu) . (7) 

N,N'=O..A 

The Fourier transform of this quantity is then fed to a suitably adapted version of WaveBurst, 



one of the burst event trigger generators used in the LIGO data analysis 19|, |20]. The 



WaveBurst algorithm make use of the wavelet decomposition, and among the bank of wavelet 
packets we picked the Symlet base with filter length sixty [24]. Using an orthogonal wavelet 
transformation the time series are converted into wavelet series Wij where i is the time 
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index and j the wavelet layer index. Each wavelet layer can be associated with a certain 
frequency band of the initial time series. The time frequency resolution varies according to 
the decomposition level, if A/" is the number of data in the initial time series, at decomposition 
level n the number of layers is 2", each with Af x 2^" points (the procedure can go on as 
long as A/'2~" is an integer). 

For each layer a fixed fraction of the coefficients W^jS with the largest absolute magni- 
tude is selected and WaveBurst keeps as triggers only those coefficients whose selection is 
robust under change of the decomposition level. Each trigger consists of highlighting a con- 
nected region in the time-frequency plane where the wavelet coefficients exceed an adaptive 
threshold (e.g. the 0.5% higher coefficients). 

An eventual GW signal impinging on the detector will have different strength in different 
channels, depending on the arrival direction and polarization, whereas in H, see eq.([6]), it 
will imprint almost exactly the same signal, no matter its polarization nor arrival direction. 

Once the trigger has been established the analysis is performed on the modes /liv's, by 
collecting the values of the wavelet coefficients for each channel and for each trigger. At this 
point we can try to reconstruct the arrival direction in (at least) two different ways. 

As a first method one hand one can find the values of 6, (j) which maximize the a posteriori 
probability of having our stretch of data given that a GW hit the detector see sec JIIIBl 
However this method will give a determination of the would be GW direction no matter 
if a real GW has excited the detector or a glitch, say, or any other noise excitation not 
compatible with a transverse, traceless GW event has taken place. To confirm this direction 
determination we combine it with a different, geometric method, see sec. IIII CI When the 
two methods do not determine the same directions, within some tolerance to be discussed 
quantitatively in sec. llVt the event can be discarded as spurious. 



The posterior probability of having a given stretch of data {hj^} assuming a GW with 
polarization strength and /ix from the direction identified by the usual polar angles 9, 
is denoted hy p{{hM}\h^,h 6, (p). The mode response to the gravity wave is 



B. Likelihood method for direction reconstruction 




(8) 
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where F^'^ is the pattern function of the mode for the +, x polarization. 

Following the standard procedure 2l|, |22| , the following likelihood ratio can be defined 

p{{hr,}\H^J 



A 



(9) 



p{{hN}\Ho) ' 

where the H^^^ is the hypothesis that a GW characterized by is present in the data and 
Hq is the hypothesis that no GW is in the data. 

For stationary, Gaussian, white noise with zero mean, and by taking into account that 
the noises in the different channels are correlated, the above mentioned probability densities 
are 



J9({MI^C^.) oc exp[-((/i-0-(/^-e)W2] 
pi{hM}\Ho) (X exp[-ih-h)s/2] , 

and the logarithm of the likelihood ratio can then be expressed as 



£ = ln(A) = 7^e[(^%]--(^05• 



(10) 



(11) 



By maximizing eq. (fTTi) with respect to h+, hx a function L{9, 0) of the angles alone is found 

(12) 



L{9, 



(F- . F-)sih . + (F+ ■ F+)sih - FX)! - 2(/. ■ F+)sih • F>^)siF+ ■ F>^)s 

The values of 6*, (which enter the expression for the F^r's ) maximizing J2{trigger}L{9,4>), 
that is the sum of the likelihood over every point of time-frequency plane exceeding the 
threshold, give the arrival direction of the candidate event. 

C. Geometric method for direction reconstruction 



Another method to reconstruct the direction of arrival of the GW by exploiting the 
(redundancy of the) five quadrupolar modes is based on linear algebra considerations [lO|| . A 
general metric perturbation hij (the temporal components are suppressed here) representing 
a GW coming from the z axis can be written as 

f h+ hy, 0^ 



h 



hy, — /i+ 

\ J 



(13) 
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For a different incoming direction the specific shape of the metric perturbation will be 
different, but the eigenvalues of the matrices will be the same, thus implying that a generic 
GW is always represented by a traceless matrix with a zero eigenvalue (thus implying that 
the non zero eigenvalues are equal in magnitude and opposite in sign). In particular the 
direction of the zero eigenvalue is the propagation direction of the GW. Of course none of the 
eigenvalues is expected to be exactly zero at every instant within the duration of the trigger, 
so for each instant of time the direction of the smallest absolute magnitude eigenvalue is 
computed. We then weight each direction with an empirical factor 1/r, taking into account 
how the eigenvalues are close to the ideal, noiseless situation described by ( fT3l) . 

After ordering the three eigenvalues Aj so that |Ao| < |Ai| < IA2I we define the following 
quantities 

V2\\o\ 



(14) 



For a perfect GW like f|T3l) r vanishes, thus the smallest it is, the less the noise is contami- 
nating the GW signal. 

IV. RESULTS 

To check our method we injected in software both GW signals and signals exciting of 
equal strength on the different channels, corresponding to a metric perturbation of the type 

^(y3-l)/y3 1 1 \ 

Kj= 1 _(73 + i)/^ 1 , (15) 

^ 1 1 

which is incompatible with a GW. We investigated how well we could recover the injections 
of both real and fake signals in each of the channels with different waveforms. 

We injected sine-gaussian with a time-width of 50msec, central frequency 2825Hz and 
with hrss = 10"2°Hz"^/^ see fig. m Fig. prefers to an injection of 25 polarized GW-signals 
incoming from the direction specified by the angle 6 = 55°, (p = 45°, whereas fig. prefers to a 
signal of the type described by eq. (fT5|l . or fake- injection, with the same shape and strength. 
The two direction identification methods agree much better in the case of GW rather than 
in the fake-injection case see fig. U\ where an almost flat distribution of distances is obtained. 
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Figure 4: The stars indicate the injections, while the dashed line is the strain sensitivity for 
the arrival direction of the injections. For comparison, we have displayed also also the averaged 
sensitivity (continuous line), as well as the ones corresponding, at any given value of the frequency, 
to the best and worst possible directions. Finally, the thick dashed line is the LIGO current 
averaged on the sky. 




Figure 5: The dots indicate direction reconstruction through the "determinant" method, the crosses 
through the likelihood one. The injection direction is marked by a circle. Both dots and crosses 
fall very well on the top of the injection direction. 
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Figure 6: The dots indicate direction reconstruction through the "determinant" method, the crosses 
through the Ukehhood one. The signal injected does not correspond to a GW. 
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Figure 7: Distribution of angular distances between the directions reconstructed with the two 
different methods in the case of GW-injections (grey-filled) and fake- injections (black-transparent). 



We then propose the combination of these methods can be used to reduce the false alarm 
rate of a single detector. 

Still, from the figH it can be seen that by setting the maximum acceptable distance to 
0.1, say, one has perfect efficiency and still picks in some 8% of the non GW-signals, thus 
suggesting that operating in coincidence with another detector is advisable. Moreover fig. [8] 



hints that the 
the direction. 



ikelihood method is more efficient than the determinant one in determining 



25| 
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Figure 8: Distribution of angular distances between the injection direction and the directions re- 
constructed with the hkehhood (grey-filled) and with the determinant method (black-transparent) 
in the case of GW-injections. The average in the two cases is respectively = 0.030 {a = 0.015), 
Hd = 0.037 {a = 0.023). 

We thus conclude these display of preliminary results by claiming that the information a 
spherical detector can provide can indeed determine the arrival direction of a GW and can 
discriminate between real GW's and glitches, say, even if more work is needed to make this 
result more quantitative and solid. 



V. CONCLUSION 



We have simulated a spherical detector and shown that a resonant sphere detector is 
capable of detecting the direction of arrival of a gravitational wave. We have shown that at 
(amplitude) SNR~ 53 (corresponding to signals with hrss = 1 x 10~^°Hz""^/^) the hkelihood 
method is well able to determine the direction, but it is not able to discriminate between 
real GW events and other kind of excitations, as for instance one with all modes excited 
equally. By cross-checking this method with a different one like we did it is possible to 
obtain an independent determination of the direction, which leaves substantially unaltered 
the efficiency and lower the false alarm rate to 8% of the rate of fake events. More work has 
to be done to assess the method at different at SNR and with higher statistics. 
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To further lower the false alarm rate it is advisable to work in coincidence, with an- 
other sphere, which could give an additional determination of the arrival direction, or with 
interferometers, which are poor in direction source determination, even if used as a network. 
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